% fig2_all.m ——— 三联图复刻论文 Fig.2 ———
clear; clc; close all;

%% —— 全局参数 —— 
beta   = 1;      % 论文中 β
xi     = 0.8;    % 论文中 ξ
mu     = 1;      % 论文中 μ
A      = 3.4;    % 外激励幅值 A
omega  = 1;      % 外激励角频率 ω

%% —— 时间控制参数 —— 
dt    = 0.01;    % 时间步长
T0    = 1000;    % 丢弃瞬态时间
Tsim  = 1200;    % 统计时间
N0    = floor(T0/dt);
Nsim  = floor(Tsim/dt);

%% —— 扫描 α 范围 —— 
alpha_min  = 1.0;
alpha_max  = 4.0;
Nalpha     = 300;
alpha_list = linspace(alpha_min, alpha_max, Nalpha);

%% —— 数据预分配 —— 
peak_x   = cell(1, Nalpha);     % Fig2(a) 峰值
LEs_list = zeros(1, Nalpha);    % Fig2(b) 最大李指数
H_avg    = zeros(1, Nalpha);    % Fig2(c) 平均哈密顿能量

%% —— 并行循环计算三组数据 —— 
parfor ia = 1 : Nalpha
    alpha = alpha_list(ia);

    % ——— (1) 仿真并提取 x 峰值 —— 
    S = [0.2;0.1;0.2];
    x_ser = zeros(N0+Nsim,1);
    for k = 1:(N0+Nsim)
        t = (k-1)*dt;
        mu_s = A*cos(omega*t);
        % RK4 一步
        k1 = f(S, alpha, beta, xi, mu, mu_s);
        k2 = f(S+0.5*dt*k1, alpha, beta, xi, mu, mu_s);
        k3 = f(S+0.5*dt*k2, alpha, beta, xi, mu, mu_s);
        k4 = f(S+    dt*k3, alpha, beta, xi, mu, mu_s);
        S = S + dt*(k1+2*k2+2*k3+k4)/6;
        x_ser(k) = S(1);
    end
    pks = findpeaks(x_ser(N0+1:end));
    peak_x{ia} = pks;

    % ——— (2) 最大李雅普诺夫指数 —— 
    Ls = LEs(alpha, beta, xi, mu, A, omega, dt, T0, Tsim);
    LEs_list(ia) = Ls;

    % ——— (3) 平均哈密顿能量 —— 
    S = [0.2;0.1;0.2];
    Hsum = 0;
    for k = 1:(N0+Nsim)
        t = (k-1)*dt;
        mu_s = A*cos(omega*t);
        % RK4
        k1 = f(S, alpha, beta, xi, mu, mu_s);
        k2 = f(S+0.5*dt*k1, alpha, beta, xi, mu, mu_s);
        k3 = f(S+0.5*dt*k2, alpha, beta, xi, mu, mu_s);
        k4 = f(S+    dt*k3, alpha, beta, xi, mu, mu_s);
        S = S + dt*(k1+2*k2+2*k3+k4)/6;
        if k > N0
            xk = S(1); yk = S(2);
            Hsum = Hsum + ( xk^2/(2*alpha) + yk^2/2 );
        end
    end
    H_avg(ia) = Hsum / Nsim;
end

%% —— 绘图 —— 
figure;

% —— (a) x 峰值分岔图 —— 
subplot(3,1,1);
hold on;
for ia = 1:Nalpha
    scatter( repmat(alpha_list(ia), size(peak_x{ia})), peak_x{ia}, 5, 'b.' );
end
xlabel('\alpha'); ylabel('x 峰值');
title('(a) x_peak 分岔图');
xlim([alpha_min alpha_max]); grid on;

% —— (b) 最大李指数谱 —— 
subplot(3,1,2);
plot(alpha_list, LEs_list, 'b-', 'LineWidth',1.2);
hold on; plot(xlim, [0 0], 'k--','LineWidth',1);
xlabel('\alpha'); ylabel('最大李雅普诺夫指数');
title('(b) LE 指数图');
xlim([alpha_min alpha_max]); grid on;

% —— (c) 平均哈密顿能量 —— 
subplot(3,1,3);
plot(alpha_list, H_avg, 'r-', 'LineWidth',1.2);
xlabel('\alpha'); ylabel('\langle H \rangle');
title('(c) 平均哈密顿能量');
xlim([alpha_min alpha_max]); grid on;
